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Abstract. A Lorenz map is a Poincare map for a three-dimensional Lorenz 
flow. We describe the theory of renormalization for Lorenz maps with a critical 
point and prove that a restriction of the renormalization operator acting on 
such maps has a hyperbolic fixed point. The proof is computer assisted and 
we include a detailed exposition on how to make rigorous estimates using a 
computer as well as the implementation of the estimates. 



1. Introduction 

The theory of renormaUzation has since its introduction about thirty years ago 
become a central tool in the study of dynamical systems. It is used, roughly speak- 
ing, to analyze maps having the property that the first-return map to some small 
part of the phase space resembles the original map itself. This property is usually 
associated with maps which lie at the "boundary of chaos", like the prototypi- 
cal example of a unimodal map at the end of a period-doubling cascade. Such 
period-doubling cascades have been observed for maps as well as flows, but most 
renormalization results so far have been for one-dimensional maps (unimodal and 
circle maps, see e.g. [5] and |17j ) with some results for higher dimensional maps 
(Henon maps, see f3^ and [1]). This paper contains new results for a class of so-called 
Lorenz flows whose dynamics can be described using the theory of renormalization. 

A Lorenz flow is a three-dimensional flow possessing a singularity of saddle type 
with a one-dimensional unstable manifold intersecting the two-dimensional stable 
manifold. If the Poincare map to a surface straddling the stable manifold can be 
foliated in such a way that the leaves are invariant and contracted exponentially 
by the Poincare map, then the dynamics of the flow is determined by the one- 
dimensional map induced by the action of the Poincare map on the leaves (see 
Chap. 14 of [3)- Such one-dimensional Lorenz maps are increasing with a jump 
discontinuity at the point corresponding to the stable manifold. They have been 
studied extensively under the additional assumption that they be expanding (see in 
particular [6] and |16j). but a much wider variety of dynamics is exhibited if there 
is also some contraction present in the form of a critical point (see |13j and [9] ) and 
this is the situation we will consider. 

The main result of this paper is that a restriction of the renormalization operator 
on the space of Lorenz maps has a hyperbolic fixed point, which is proved using 
the contraction mapping theorem on an associated operator. We use a computer 
to rigorously compute estimates that shows that this associated operator is indeed 
a contraction. This method was pioneered by Lanford fTT| (see also [T^) when he 
proved the existence of a fixed point of the period-doubling operator on unimodal 
maps. However, Lanford's paper only gives a brief outline of the method he employs 
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without an actual proof, so we have gone through quite a lot of pains to include all 
the missing details in this paper (many of which were borrowed from |10j). 

This article is roughly divided into two halves: In Section [5] we give all the 
necessary definitions to state the renormalization conjecture and the main theorem, 
and then we go on to prove several consequences of the main theorem in Section |3l 
In Section|4]we describe the method used to prove the main theorem and in Section[5] 
we give the proof. This concludes the first half. 

In the second half of the paper we give exact details on how the estimates needed 
to prove the main theorem are implemented on a computer^ The literature on this 
type of computer assisted proof seems to have a tradition of never including these 
details, most likely because it would require an order of thousands of lines of source 
code. We make a conscious break from this tradition and show how to implement 
all estimates in only 166 lines of source codeQ The key behind this reduction in 
size is to use a functional programming language since it allows us to program in 
a declarative style: we specify what the program does, not how it is accomplished. 
This also has the benefit that functions cannot have side- effects (the output from 
a function only depends on its input) which makes it easier to reason about the 
source code. In our context this is extremely important since it means that we 
can check the correctness of each function in complete isolation from the rest of the 
source code (and a typical function is only one or two lines long which simplifies the 
verification of individual functions). To further minimize the risk of programming 
errors we choose a strongly typed language since these are good at catching common 
programming errors during compilation. 

We take this opportunity to advocate the programming language Haskell for 
tasks similar to the one at hand — it has all the benefits mentioned above and 
more, but at the same time manages to produce code which runs very fast (thanks 
to the GHC compiler). Unfortunately, many readers will probably have had little 
prior exposure to Haskell and for this reason we have in Appendix |E] included a 
brief overview of Haskell as well as a table highlighting its syntax to aid the reader 
in understanding the source code. 

2. Statement of the main result 

In this section we state the main result, but in order to do so we first need quite 
a few definitions. 

Definition 2.1. A Lorenz map / on a closed interval I = [l,r] {I < < r) is a 
monotone increasing continuous function from / \ {0} to / such that /(O^) — r, 
/(0+) = I (i.e. / has a jump discontinuity at 0)0 

We require that f{x) = ip{\x\'') for aU x G {I, 0), where tpis a symmetric^ analytic 
map defined on some complex neighborhood of [1,0], and similarly f{x) = 'ip{\x\'') 
for x € (0,r), where ip is a symmetric analytic map defined on some complex 

^ We have structured the article so that a reader with no interest in the details of the computer 
estimates can skip the second half of the article. 

^ This includes: definition of main operator and its derivative (40 lines), an interval arithmetic 
library (30 lines), a library for computing with analytic functions (65 lines), a linear equation 
solver (15 lines). 

The notation /(0~) is shorthand for lim^j-i-o f{x) and this limit is assumed to exist; /(O"*") is 
defined analogously as the right-hand limit. 
^ Here 'symmetric' means ip{z) = <^(z). 
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Figure 1. A Lorenz map (the fixed point of Tlieorcm l2.1ip . 

neighborliood of [0,r]. Tlie maps tp and -0 are called the analytic parts of /. The 
constant p > 1 is called the critical exponent of / (and is independent of /). 

Assume / is defined on [—1, r] and let 3 be a Lorenz map on [—1, r'] with analytic 
parts ((/?', V')- We define a metric on the set of Lorenz maps by 

d{f,g)=max{\\^-^'\\,U-tP'\\}, 

where ||-|| denotes the usual sup~norm on analytic functions. (For Lorenz maps 
with different domains we first perform a linear coordinate change to ensure that 
their domains are of the above form, then apply the above formula for the metric.) 

Remark 2.2. The condition p > 1 ensures that Df{x) ^ as a; — > (from the 
left or the right) and for this reason we call the critical point. Because of the 
discontinuity at there are two critical values, namely I and r. 

The smoothness required in our definition of Lorenz maps is not essential for a 
satisfactory renormalization theory, but our current results are only in this category 
(which is not a big restriction since they can most likely be extended to C for r ^ 3 
along the lines of 4 and [5]). For a discussion on what to expect when the minimum 
smoothness threshold is approached from below, see |2]. 

Definition 2.3. A branch of /" is a maximal open interval B on which /" is 
monotone (here maximality means that if A is an open interval which properly 
contains B, then /" is not monotone on A). 

To each branch B of /" we associate a word w(i?) = {ao, . . . , an-i} on two 
symbols by 

^ ^ fo ifP(S)c(/,0), 
\l ifP(S)c(0,r), 

for j = 0, . . . , n — 1. 

Definition 2.4. A Lorenz map / on / is renormalizable if there exists a maximal 
interval U <^ I containing such that the first-return map / to C/ is a Lorenz map 
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on U. In this situation we define the renormalization TZf of / as the first-return 
map rescaled via an increasing hnear map h : I ^ U which takes to itself and 
the left endpoint of / to the left endpoint of [70 

TZf = h-^ofoh. 

The operator TZ is called the renormalization operator. 

Remark 2.5. When defined on the space of Lorenz maps with analytic branches 
the renormalization operator is differentiable and its derivative is a compact linear 
operator. This follows from the fact that TZf only evaluates / on a strict subset 
of the domain of / (see Sections 19.41 and 19. 5|) . On the other hand, if we only were 
to demand C-smoothness for the branches of our Lorenz maps then TZ would no 
longer be differentiable (see [S] and [TH Ch. VI. 1.1]). 

Definition 2.6. Let / be a renormalizable Lorenz map with associated first-return 
map / and return interval U. Then there exists (unique) integers a,b > 2 such that 
/li = P and f\R = /^ where L = U n (1,0) and R = U n (0,r). The interval 
L is contained in a branch A of /° with associated word a = u}{A), and similarly 
i? C i? for a branch B of with /3 — uj{B). The pair of words {a, [3) is called the 
type of renormalization. 

The notation TZa,i3 will be used to denote the restriction of TZ to the set of Lorenz 
maps which have renormalizations of type (a,/3). 

Remark 2.7. In the kneading theory for unimodal maps a finite part of the knead- 
ing sequence of a point does not contain enough information to recover the exact 
ordering of the corresponding points on the real line. This leads to the introduction 
of "unimodal permutations" to describe the combinatorics of unimodal renormal- 
ization. Since Lorenz maps are strictly increasing and have no fixed points the 
ordering of points of a finite part of an orbit is completely determined by the cor- 
responding kneading information so there is no need for this extra complication. 
However, we may still ask which pairs of words (a, /3) give rise to valid types of 
renormalization. The answer to this question is given in T5|; we will not go into 
details as this would require more definitions, but suffice to say that there is a 
simple admissibility condition stated in terms of the shift operator acting on words 
on two symbols. 

Definition 2.8. Let / be a Lorenz map such that TZ"f is defined for every pos- 
itive integer n. In this situation we say that / is inGnitely renormalizable. The 
combinatorial type of / is the sequence of words {(ao, /3o), («!, /^i), • • • }, where / 
has renormalization of type (ao, /3o), T^f has renormalization of type (ai, /3i), and 
so on3 If the length of the words ak and are bounded in k then / is said to be 
of bounded combinatorial type. 

Remark 2.9. Let / be of combinatorial type {ctq, cri, . . . }, where ak = (afc,/3fc). 
Then TZf has combinatorial type {cti, CT2, . . . }. In other words, TZ shifts the combi- 
natorial type to the left. 

^ The choice of rescaling is somewhat arbitrary (as long as it is affine) — we have chosen it so 
that the critical point and the left endpoint of the domain of / are fixed under renormalization, 
whereas the right endpoint may move. Another natural choice is to fix the endpoints of the domain 
of / but then the critical point may move. 

^ When talking about the combinatorial type we implicitly assume that it is admissible in line 
with the discussion of Remark l2.7l 
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We are now ready to state the main result, but before doing so we mention 
the renormahzation conjecture for Lorenz maps (the statement is taken from the 
corresponding resuh for unimodal maps, see [5]): 

Conjecture 2.10 (Renormahzation horseshoe). The limit set of TZ acting on the 
space of Lorenz maps of hounded combinatorial type is a hyperbolic Cantor set where 
TZ is conjugate to the full shift in a finite number of symbols, for every critical 
exponent p > I. 

The above theorem represents an ultimate goal for the theory of renormahzation. 
However, our results are much more modest in that we only prove that locally at 
one point in the limit set of TZ the above conjecture holds: 

Theorem 2.11 (Main Theorem). Let a = {0, 1} and P ^ {1,0,0}. The restricted 
renormahzation operator TZa.p acting on the space of Lorenz maps with critical 
exponent p — 2 has a hyperbolic fixed point. 

Proof. This is a direct consequence of the estimates in Theorem 15.31 and the dis- 
cussion in Section 14.11 □ 

We would like to address the question as to how far towards the renormahza- 
tion conjecture our method of proof can take us. Unfortunately the answer is "not 
very" . Theoretically, given a critical exponent p > 1 and any periodic combinato- 
rial type {ao, . . . , cr„, co: ■ • ■ : cr„, . . . } we could check estimates similar to those of 
Theorem l5.3l in order to deduce the existence of a hyperbolic fixed point. However, 
implementing these checks even for the simple combinatorial type at hand requires 
a significant effort so this does not really have any practical significance. Despite 
these shortcomings we still think that our current result is an important first step 
in the theory of renormahzation of Lorenz maps. 

It is also interesting to ask if any of the methods from the theory of renormal- 
ization of unimodal maps can be used to prove the renormahzation conjecture for 
Lorenz maps. We do not know the answer to this question but it seems unlikely 
since the unimodal theory is based on complex analytic methods that do not have 
any obvious generalization to Lorenz maps (since Lorenz maps have a point of 
discontinuity). For this reason the renormahzation theory for Lorenz maps poses 
new and significant difficulties. However, we can use some results from unimodal 
renormahzation as the following remark shows. 

Remark 2.12. The fixed point of Theorem 12. lll is the simplest non-unimodal fixed 
point of TZ. By this we mean that if a and (3 both have length 2 (i.e. a — {0, 1}, 
/? — {1, 0}) then the fixed point of the period-doubling operator on unimodal maps 
corresponds to a fixed point for TZa^i3 as follows. 

Let g : [—1,1] — )• [—1,1] be the fixed point for the period-doubling operator 
normalized so that g{0) = 1. Then g is an even map that satisfies the Cvitanovic- 
Feigenbaum functional equation 



Now define a Lorenz map / by /|[-i,o) ~ 9 and /|(o,i] — — 3- It is easy to check 
that the first-return map f to U ~ [—A, A] is / = and that U is maximal. Thus 



X~'g'{Xx), 



A = -5(1)- 
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Figure 2. Illustration of the dynamical intervals of generations 
0, 1, 2 for renormalization of type a = {0, 1}, /? = {1, 0, 0}. 



which shows that / is a fixed point of TZa,p- 

3. Consequences of the main result 

The existence of a hyperbolic renormalization fixed point has very strong dynam- 
ical consequences, some of which we will give a brief overview of here. Through- 
out this section let TZ denote the restricted renormalization operator TZa, where 
a = {0, 1} and (5 = {1,0,0}, and let denote the fixed point of Theorem 12. 11 1 

Corollary 3.1 (Stable manifold). There exists a local stable manifold Wf^j, at 
consisting 0/ maps in a neighbourhood of which under iteration of TZ remain in 
this neighbourhood and converge with an exponential rate to 

The local stable manifold extends to a global stable manifold consisting of 
maps which converge to under iteration of TZ. If f £ W"* then f is infinitely 
renormalizable. 

Proof. The existence of a stable manifold is a direct consequence of the stable and 
unstable manifold theorem. If / converges to then 7?."/ is defined for all fc > 0, 
which is the same as saying that / is infinitely renormalizable. □ 

We now turn to studying the dynamical properties of maps on the stable mani- 
fold. Let / e ^ then the times of closest return (a„, 6„) are given by the recursion 

Oji+l = On + bn, 0,1 = 2, 

bn+i = 2a„ + 6„, bi = 3, 

These determine the first-return interval [/„ = cl{L„ U Rn} for the n-th renormal- 
ization, by 

in = (/''"(o+),o), i?„ = (o,r"(o-)). 

In other words: the first-return map /„ to Un is given by fn{x) — {x) \i x ^ Ln 
and fn{x) = f''"(x) if x e i?„. 
Define 

Li^f^{Ln), fc = 0,...,a„-l, 
R^ = f''iRn), fc = 0,...,&„-l. 

The collection of these intervals (over fc) form a pairwise disjoint collection for 
each n, called the intervals of generation n (see Figure [2|). 

Theorem 3.2 (Cantor attractor). /// G W then the closure of the critical orbits 
of f is a measure zero Cantor set Aj which attracts almost every point in the 
domain of f . 



A RENORMALIZATION FIXED POINT FOR LORENZ MAPS 



7 



Proof. The critical orbits form the endpoints of the dynamical intervals {L'^} U 
{Rn}, so A/ is contained in 

(1) fl c\{E„ UF„}, where E„=[jL'^ and F„ = |J i?^ 

" k k 

Note that En+i C £'„ and i^„+i C F„. 

First assume that / = /*. Then the first-return maps /„ are all equal to / itself 
(up to a linear change of coordinates) so the total lengths of iJ„ and F„ shrink 
with an exponential rate (the position of f/„+i inside Un is the same for all n, so 
we can apply the Macroscopic and Infinitesimal Koebe principles as in the proof of 
the "real bounds" in [Ml Ch. VI.2]). Hence the intersection ([T]) is a measure zero 
Cantor set, and consequently A/ is as well. 

Now, if / is an arbitrary map in then TZf converges to fi,. In other words, 
the first-return maps {/«} converge to (up to a linear change of coordinates). 
Now use the same arguments as above. 

For a proof that A/ is an attractor, see [15]. □ 

Theorem 3.3 (Rigidity). If f,g G W"* then there exists a homeomorphism h : 
A/ — S> Ag conjugating f and g on their respective Cantor attractors. If furthermore 
f,g £ y^ioc' then h extends to a C^"*"" dijfeomorphism on the entire domain of f . 

Proof Define /i(/"(0-)) = g"(0-) and /i(/"(0+)) = g"(0+). This extends contin- 
uously to a map on A/ as in the proof of [Ml Proposition VI. 1.4]. 

If /,5 G Wf„^ then there exists C > and A < 1 such that d(/" ,5") < CA" so 
we can use an argument similar to that in |14[ Theorem VI. 9. 4] to prove the second 
statement. □ 

Remark 3.4 (Universality). The second conclusion of Theorem 13.31 is a strong ver- 
sion of what is known as "metric universality" : the small scale geometric structure 
of the Cantor attractor does not depend on the map itself (only on the combinato- 
rial type and the critical exponent). That is, if we take two maps /, 5 G Wj^^, and 
zoom in around the same spot on their Cantor attractors then their structures are 
almost identical since a differentiable map (i.e. the extended h) is almost linear if 
one zooms in closely enough. 

For example, the limit of [L„_|_i|/|L„| as n — >■ 00 exists and is independent of / (it 
equals the ratio |L2(/*)l/|-^i(/*)| for /*). More generally, the multifractal spectrum 
(and Hausdorff measure in particular) of A/ does not depend on / (only on fi,). 

We also want to mention another type of universality called "universality in the 
parameter plane" where the unstable eigenvalues of DTf^ governs the structure of 
the parameter plane for families of Lorenz maps. However, in order to present the 
details we would need more definite information on the structure of the spectrum 
of DTf^ so we will have to return to this discussion in another paper. 

4. Outline of the computer assisted proof 

In this section we give a brief outline of the method of proof and how to calculate 
rigorous estimates with a computer. 

4.1. Method of proof. Given a Frechet differentiable operator T with compact 
derivative on a Banach space X of analytic functions we would like to prove that 
T has a hyperbolic fixed point. The main tool is the following consequence of the 
contraction mapping theorem: 
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Proposition 4.1. Let ^ be a Frechet differ entiable operator on a Banach space X , 
let fo G X, and let Br{fo) C X be the closed ball of radius r centered on /q. // 
there are positive numbers e, 9 such that 

(1) \\D^f\\<9,forallf€Br{fo), 

(2) ||<f/o-/o|| <e, 

(3) £ < (1 - ey, 

then there exists f-i, € B^ifo) such that $/^, = f^ and $ has no other fixed points 
inside i?r(/o)- Furthermore, ||/* — /o|| < e/{l — 0). 

Our strategy is to find a good approximation fo of a fixed point of T and then use 
a computer to verify that the conditions on r,e,9 hold if r is chosen small enough. 
Unfortunately, this is not possible for T itself since in our case it is not a contraction 
(it has expanding eigenvalues) so first we have to turn T into a contraction without 
changing the set of fixed points. This is done by using Newton's method to solve 
the equation Tf — / = 0, which results in the iteration 

f^f~{DTf-I)-\Tf^f), 

where / denotes the identity operator on X. The operator we use is a slight sim- 
plification of this, namely 

*/ = / - (r - i)-\Tf - /) = (r - i)-\r - T)f, 

where F is a finite-rank linear approximation of DTf^ (chosen so that F — / is 
invertible). The operator $ is a contraction if /o and F are chosen carefullyQ 
Note that $/ = / if and only if Tf = /, so once we verify that the conditions of 
Proposition 14. 1 1 hold for $ it follows that T has a fixed point. 

To prove hyperbolicity we need to do some extra work. The derivative of $ is 

D<^f^{r-i)-\r-DTf). 

At this stage we will have already checked that the norm of this is bounded from 
above by 1. By strengthening this estimate to 

||(F - e'*/)-i(F - DTf)\\ < 1, yte M, V/ e i3.(/o), 

we also get that DTf^ is hyperbolic at the fixed point fi,. To see this, assume that 
e'* is an eigenvalue of DTf^ with eigenvector h normalized so that \\h\\ = 1. Then 

\\(r-e^'l)-\T-DTjJh\\ = ||(F-e^*/)-i(F-e^*/)/.|| = ||;.|| ^ 1, 

which is impossible. Since DT was assumed to be compact we know that the spec- 
trum is discrete, so the lack of eigenvalues on the unit circle implies hyperbolicity. 

4.2. Rigorous computer estimates. In order to verify the above estimates on a 
computer we are faced with two fundamental problems: (i) arithmetic operations on 
real numbers are carried out with finite precision which leads to rounding problems, 
(ii) the space of analytic functions is infinite dimensional so any representation of 
an analytic function needs to be truncated. 

The general idea to deal with these problems is to compute with sets which 
are guaranteed to contain the exact result instead of computing with points: real 



^ Here is how to choose /q and F: use the Newton iteration on polynomials of some fixed 
degree to determine /o and set F = DTf^. The hardest part is finding an initial guess such that 
the iteration converges. 
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numbers are replaced with intervals, analytic functions are replaced with rectangle 
sets Aq X ■ ■ ■ X Ak X {C} in M" representing all functions of the form 

{aoH h afcz'' + | £ Aj, j = 0, . . . , fc, \\h\\ < C}, 

where {Aj} are intervals. This takes care of the truncation problem and the round- 
ing problem is taken care of roughly by "rounding outwards" (lower bounds are 
rounded down, upper bounds are rounded up). Once these set representations have 
been chosen we lift operations on points to operations on sets. Since the form of 
these sets are most likely not preserved by such operations, this lifting involves 
finding bounds by sets of the chosen form (e.g. if F and G are rectangle sets of 
analytic functions and we want to lift composition of functions, then we have to 
find a rectangle set which contains the set {f o g \ f E F, g G G}.) 

Section [7] contains all the details for computing with intervals and Section [9] 
contains all the details for computing with rectangle sets of analytic functions. 

Let us make one final remark concerning the evaluation of the operator norm 
of a linear operator L on the space of analytic functions. In order to get good 
enough bounds on the estimate of the operator norm we will use the ^^-norm on 
the Taylor coefficients of analytic functions. The reason for this is that estimating 
the operator norm with 

ll^ll = sup IIT/II 
ll/lt<i 

will usually result in really bad estimates. With the ^^-norm, if we think of L as 
an infinite matrix (in the basis {z*^}), the operator norm is found by taking the 
supremum over the norms of the columns of this matrix, that is 

||L||=sup||LCfe||, ^k{z)^z\ 

k>0 

Evaluating the norms of columns gives much better estimates and for this reason 
we choose this norm. See Section [9. Ill for the specifics. 



5. Proof of the main theorem 

First we restate the definition of the restricted renormalization operator, then 
we change coordinates and restate the main theorem. 

5.1. Definition of the operator. From now on we fix the domain of our Lorenz 
maps to some interval [— 1, r]. The right endpoint cannot be fixed since it generally 
changes under renormalization (we will soon change coordinates so that the domain 
is fixed). 

Instead of dealing with functions with a discontinuity we represent a Lorenz 
map F by a pair (/, g), with / : [—1, 0] [—1, r], /(O) — r, and g : [0, r] — > [—1, r], 
5(0) = -1.^ 

With this notation, the first-return map to some interval U will be of the type 
{F"-,F^)\u if F is renormalizable. For the type a = {0, 1}, /3 = {1, 0, 0}, we can 
be more precise: in this case a = 2, 6 = 3 and the first-return map is of the form 
{g ° f 1 f ° f ° g)\u if it is renormalizable. 
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Let T denote the restricted renormalization operator TZa,i3, and fix the critical 
exponent p = 2. If T{f,g) — {f,g) then T is defined by 

/(z) = A-i.g o /(Az), 

g{z) = \-'fofog{Xz), 

A = -/'(-!). 

5.2. Changing coordinates. To ensure the correct normalization (5(0) = —1) 
and the correct critical exponent {p = 2) we make two coordinate changes and cal- 
culate how the operator T transforms. We will also carefully choose the domain of 
T so that all compositions are well-defined (e.g. Xz is in the domain of / etc.). This 
is checked automatically by the computer (and also shows that T is differentiable 
with compact derivative, since / and g are analytic). Finally, it is important to 
realize that the choice of coordinates may greatly affect the operator norm of the 
derivative; not every choice will give a good enough estimate. 

The domain of T is chosen to be contained in the set of Lorenz maps (/, g) 
with representation f{z) — (j){z'^) and g{z) — 'ip{z'^), where ip and 'ip have domains 
{z : \z — 1\ < s} and {z : \z\ < t}, respectively (the constants s and t will soon be 
specified) . Rewriting T in terms of (j) and "0 gives 

0(z) = A-V(0(A2^)2) 

i,{z) = \-^m^{^^^ff) 

A = -0(0(1)2) 

This coordinate change ensures the correct critical exponent. 

The next coordinate change is to fix the normalization and also to bring the 
domain of all functions to the unit disk. Fixing the normalization has the benefit 
that the error involved in the evaluation of A is minimized (since we only need to 
evaluate / close to z — see Section [978]) . Changing all domains to the unit disk 
simplifies the implementation of the computer estimates. 

Definition 5.1. Define X to be the Banach space of symmetric (with respect to 
the real axis) analytic maps on the unit disk with finite £^-norm. That is, if / £ X 
then /(z) = ^flfcz'^ with e M and ||/|| — X^l^fcl < oo- 

Definition 5.2. Define Y = XxX with the norm ||(/,5)l|y = + hWx and 

with linear structure defined by a{f,g) + (3{f',g') = {af + f3f', ag + (3g')- Clearly 
y is a Banach space (since X is). 

Change coordinates from 0, ip to (/, g) ^ Y (note that / and g are not the same 
as above) as follows 

0(z)=/([z-l]/s), 
^{z) = -1 + z-giz/t), 
where we will choose s = 2.2 and t — 0.5. Rewriting T in terms of / and g gives 

/(-) = {-1 + ./ (A^ h + - i)' • (i • / (A^ h + i] - IT) } , 

9H = ^ {1 + A- V (1 . / {Xhwg (X^w) ■ [Xhn^g (X^u.) - 2] • 1)^ - l) } , 

\ = -.f{lf{Of-l]/s). 
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This is the final form of the operator that will be studied. 

5.3. Computing the derivative. In order to simplify the computation of the 
derivative of T we break the computation of T down into several steps as follows: 



Pf{w) 


^ \^ ■ {w + s-^) - s-^ 


Pg{w) 


= X^w 


h 


= f°pf 


91 


= 9°Pg 


h 




92 


^t-Pg ■ 91 


h 


= hit 


93 


= 92 ■ (.92 - 2)/s 


h 


= 9° fs 


9i 


= / °53 


h 




95 


= (54 - 


k 




9e 


= / 055 






97 


= 9e/X 






58 (w) 


= igriw) + l)/(t • w) 



With this notation we have that T(/, g) = (/e, 58)- Note that the result of 37 (w) + 1 
is a function with zero as constant coefficient so in the implementation of gs we will 
not actually divide by w, instead we will 'shift' the coefficients to the left. 

It is now fairly easy to derive expressions for the derivative. If / is perturbed by 
Sf and g is perturbed by Sg, then the above functions are perturbed as follows: 



Spf{w) ^2- X-SX-{w + s-') Spg{w) = 2-X-6X- 



w 



6fi ^ Df opf ■ dpf + df o pf 5gi ^ Dgopg ■ Spg + Sgopg 

5f2 = 2/1 (5 /i (5.92 = t ■ {6pg ■ .91 + pg ■ Sgi) 

•5/3 = Sf2/t Sgs = Sg2 ■ (52 - 2)/s + 32 • Sg2/ s 

(5/4 = Dg o /g • Sfs + Sg o /g 5g4 = Df og^- Sg^ + Sf o 53 

(5/5 = Sf2 ■ fi + h- 5fi 3g5 = 2- g4- Sg^/s 

(5/6 = '5/5/A - /s • (5A/ Sge = Df og^- Sg^ + Sf og^ 

(557 = '556/A - 56 • (5A/A^ 
Sgs{w) = Sg7{w)/{t ■ w) 
With this notation we have that DTi^j g){5f ,6g) = {Sfg^Sgs). 

5.4. New statement of the main theorem. We now state the main theorem in 
the form it will be proved. The discussion in Section [4.11 shows how this result can 
be used to deduce Theorem l2.11l 

Theorem 5.3. There exists a Lorenz map Fq and a matrix T such that the sim- 
plified Newton operator ^ — (T — /)~^(r — T) lis well-defined and satisfies: 

(1) \\D<^f\\ < 0.2, for all \\F - Fo\\ < 10^^ 

(2) ||$Fo-Fo|| <5-10-9. 

(3) ||(r - e''*/)-i(r - DTf)\\ < 0.9, for all t e M, ||F - Fo|| < IQ-^. 

Proof. The remainder of this article is dedicated to rigorously checking the first two 
estimates with a computer. The third estimate is verified by covering the unit circle 
with small rectangles and using the same techniques as in the first two estimates 
to get rigorous upper bounds on the operator norm. However, we have left out the 
source code for this estimate to keep the page count down and also because the 
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running time of the program went from a few seconds to several hours (we had to 
cover the circle with 50000 rectangles in order for the estimate to work). □ 

Remark 5.4. The approximate fixed point Fq and approximate derivative T at the 
fixed point are found by performing a Newton iteration eight times on an initial 
guess (which was found by trial-and-error) . We will not spend too much time 
talking about these approximations but they could potentially be used to compute 
e.g. the HausdorfF dimension of the Cantor attractor of maps on the local stable 
manifold. 

We did however compute the eigenvalues of F and it turns out that F has two 
simple expanding eigenvalues As ~ 23.36530 and A^, w 12.11202, and the rest of 
the spectrum is strictly contained in the unit disk. Since F is a good approximation 
of DTf^ and both operators are compact it seems clear that the spectrum of DTf^ 
also must have exactly two unstable eigenvalues. 

Lanford [12] claims that in the case of the period-doubling operator if an analog 
of the third estimate of Theorem 15.31 holds and "if F has spectrum inside the unit 
disk except for a single simple expanding eigenvalue, then the same will be true for 
DTf^." It seems plausible that a similar statement holds in the present situation 
with two simple expanding eigenvalues but have not yet managed to prove this (it 
is easy to see that if F and DTf^ were both diagonal then the third estimate would 
imply that they have the same number of unstable eigenvalues). 

6. Implementation of computer estimates 

In this section we implement the main operator and compute the estimates of 
Theorem l5.3l Before reading this section it may be a good idea to take a quick glance 
at the beginning of Section [9] in order to understand the way analytic functions are 
represented. It may also be helpful to use Table [T] in Appendix |E] to look up 
unfamiliar syntax in the source code. 

6.1. The main program. To begin with we import two functions from the stan- 
dard library that will be needed later: 

1 import Data. List (maximumBy .transpose) 

The entry point of the program is the function main, all that is done here is to 
print the result of the computations to follow: 

2 main = do putStrLn $ "radius = " ++ show beta 

putStrLn $ "|Phi(f)-f| < " ++ show eps 
putStrLn $ "iDPhil < " ++ show theta 

The initial gues^ is first improved by iterating a polynomial approximatior0 of the 
operator $ eight times (the derivative is recomputed in each iteration, so this is a 
Newton iteration): 

5 approxFixedPt = iterate (\t -> approx $ opPhi (gamma t) t) guess ! ! 8 
Compute the approximation F of the derivative DT at the approximate fixed point: 

6 approxDeriv = gamma approxFixedPt 



^See Appendix [Cl 
^See Section [9A01 
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Compute an upper bound on the distanc^I^ between the approximate fixed point 
and its image under $: 

eps = upper $ dist approxFixedPt (opPhi approxDeriv approxFixedPt ) 

Construct a balO of radius j3 centered around the approximate fixed point and 
then compute the supremum of the operator nor nO of the derivative on this balh 

theta = opnorm $ opDPhi approxDeriv (ball beta approxFixedPt) 

The rest of this section will detail the implementation of the operator $ and 
its derivative. The generic routines for rigorous computation with floating point 
numbers and analytic functions are discussed in the sections that follow. All input 
to the program (d, sf , sg, guess, beta) is collected in Appendix [Cl Instructions 
on how to run the program and the output it produces is given in Appendix [P] 

6.2. The main operator. The operator T is computed in a function called mainOp 
which takes a Lorenz map (/, g) and a sequence of tangent vectors {{Sfk,Sgk) G 
yjfe^i and returns {T{f, g), {DT(f^g)(Sfk,5gk)}k=i)- We perform both computa- 
tions in one function since the derivative uses a lot of intermediate results from the 
computation of T{f, g). 

Given a Lorenz map (f ,g) and a list of tangent vectors ds, first compute /g and 
(78 as in Section 15.31 and split the result so that the polynomial parts have degree 
at most d — \. Then compute the derivatives and return the result of these two 
computations in a pair: 

mainOp (f,g) ds = ((split d f6, split d g8) , mainDer ds) where 



1 




lambda f 








pf 




F [(l-2-l)/sf ,1'2] 


gl 




compose g pg 


pg 




F [0,1"2] 


g2 




pg * gl .* sg 


f 1 




compose f pf 


g3 




g2 * (g2 - 2) ./ sf 


f2 




f 1-2 


g4 




compose f g3 


f3 




f2 ./ sg 


g5 




(g4-2 - 1) ./ sf 


f4 




compose g f3 


g6 




compose f g5 


f5 




-1 + f2*f4 


g7 




g6 ./ 1 


f6 




f5 ./ 1 


g8 




Ishift gl ./ sg 



The actual computation of the derivative is performed next inside a local function 
to mainOp. If there are no tangent vectors, no computation is performed: 

mainDer [] = [] 

Otherwise, recurse over the list of tangent vectors and compute (J/g and 5gs, and 
again split the result so that the polynomial parts have degree at most d — 1: 

mainDer ((df,dg):ds) = (split d df6, split d dg8) : mainDer ds where 



dl 


= dlambda f df 






dpf 


= F ([2*l*dl]*[l/sf ,1]) 


dgl 


= dcompose g pg dg dpg 


dpg 


= F [0,2*l*dl] 


dg2 


= (dpg*gl + pg*dgl) .* sg 


df 1 


= dcompose f pf df dpf 


dg3 


= 2*(dg2*g2 - dg2) ./ sf 


df2 


= 2*fl*dfl 


dg4 


= dcompose f g3 df dg3 


df3 


= df2 ./ sg 


dg5 


= 2*g4*dg4 ./ sf 



lOSee Section 1931 
l^See Section 19021 
^^See Section |91T1 
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df4 = dcompose g f3 dg df3 ; dg6 = dcompose f g5 df dg5 
df5 = df2*f4 + f2*df4 ; dg7 = dg6 . /I - g6.*(dl/l~2) 

df6 = df5./l - f5.*(dl/l-2) ; dgS = Ishift dg7 ./ sg 

Note that the constants s and t of Section [531 are caUed sf and sg respectively in 
the source code. 

The above function can be used to compute the action of T by passing an empty 
hst of tangent vectors and extracting the first element of the returned pair: 

opT fg = fst $ mainOp fg [] 

Similarly, we can evaluate DT by extracting the second element: 
opDT fg ds = snd $ mainOp fg ds 

Using this function we compute an approximation F of ZJTj^ g) by evaluating the 
derivative at the 2d first basis vectorf0 of Y and approximating the result with 
polynomials and packing them into a 2d x 2d matrix (transposing the resulting 
matrix is necessary because the linear algebra routineS we use require the matrix 
to be stored in row- major order): 

gamma fg = transpose $ map (interleavePoly . approx) 
$ opDT fg (take (2*d) basis) 

Finally, the operator $ (and its derivative) is implemented by taking a Newton 
stepB with T (for convenience we pass the approximate derivative as the parameter 
m): 

opPhi m X = newton m (opT x) x 

opDPhi m X ds = [ newton m a b I (a,b) <- zip (opDT x ds) ds ] 

6.3. The rescaling factor. With our choice of coordinates the rescaling factor A 
only depends on / (and not on g): 

Xif) = -f{[f{Of-l]/s) 

The implementation is straightforward: 

lambda f = -eval f (((eval f 0)-2-l)/sf) 

If / G X is perturbed hy Sf E X then A is perturbed by SX, where 

SX^-2- ,s-i • /(O) • Sf{0) ■ Df{[J{Qf - lys) - 5f{[f{Qf - 

Derivative evaluation has to be handled carefully since we are using the ^^-norm, 
see Section |9^ If y = [/(O)^ — then y lies in the closed disk of radius |y| but 
since we need to evaluate the derivative on an open disk we first enlarge the bound 
on \y\ to get the radius /i and then evaluate the derivative on this slightly larger 
disk: 

dlambda f df = -2/sf * fO * eval df * eval (deriv mu f) y - eval df y 
where fO = eval f 

y = (f0~2 - l)/sf 
mu = enlarge $ abs y 



13 See Section [Qin 
1^ See Appendix 1X1 
15 See Section [933] 
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7. Computation with floating point numbers 

We discuss how to control rounding and avoid overflow and underflow when 
computing with floating point numbers. We show how to lift operations on floating 
point numbers to intervals and then how to bound these operations. 

7.1. Safe numbers. In order to avoid overflow and underflow during the course 
of the proof we restrict all computations to the set of safe numbers (see [10]) which 
we define as the subset of double precision floating point numbers (referred to as 
Boats from now on) x such that x — or 2~^°° < < 2^™. We say that y is a 
safe upper bound on x iS x < y (strict inequality) and y is a safe number; safe 
lower bounds are defined analogously. If a; = 0, then y — is both a safe upper and 
lower bound and there are no other safe bounds on x (this will make sense after 
reading the assumption below). 

Safe numbers allow us to perform rigorous computations on any computer con- 
forming to the IEEE 754 standard since such a computer must satisfy the following 
assumption!^ 

Assumption. Let x be a float resulting from an arithmetic operation on safe num- 
bers performed by the computer and let x be the exact result of the same operation. 
If X then either x ~ x^ or x — x~^ , where x~ is the largest float such that 
x~ < X and x^ is the smallest float such that x < x^ . Furthermore, x = if and 
only if X = 0. 

Under this assumption we know that the exact result must lie within any safe 
upper and lower bounds on x, and we know that when the computer returns a 
result of then the computation must be exact. 

Given a float x we now show how to find safe upper and lower bounds on x. 

Check if a number is safe: 

isSafe x = let ax = abs x in x == I I (ax > 2~'(-500) && ax < 2"500) 

Use this function to assert that a number is safe, abort the program otherwise: 

assertSafe x I isSafe x = x 

I otherwise = error "assertSafe: not a safe number" 

Given a float we can 'step' to an adjacent float as follows: 

stepFloat n = 

stepFloat n x = let (s,e) = decodeFloat x in encodeFloat (s+n) e 

That is, StepFloat 1 x is the smallest float y larger than x and stepFloat (-1) x 
is the largest float smaller than x, unless a; = in which case x is returned. (The 
function decodeFloat converts a float to the form s ■ 2"^, where s,e G Z, and 
encodeFloat converts back to a float.) 

Now flnding a safe upper or lower bound is easy, just step to the next float and 
assert that it is safe: 

saf eUpperBound = assertSafe . stepFloat 1 
saf eLowerBound = assertSafe . stepFloat (-1) 



^*'This statement follows from: (1) the fact that IEEE 754 guarantees correct rounding, (2) 
the result of an arithmetic operation on safe numbers is a normalized float so silent underflow to 
zero cannot occur. 
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7.2. The Scalar data type. The Scalar data type represents safe lower and 
upper bounds on a number: 

data Scalar = S [Double [Double deriving (Show,Eq) 

The first number is the lower bound, the second the upper bound. The following 
function returns the upper bound: 

upper (S _ u) = u 

We bound operations on real numbers by first lifting them to operations on 
Scalar values and then bound the resulting operations by enlarging the bound 
to safe lower and upper bounds. An operation is exact if it does not involve any 
rounding (in which case there is no need to enlarge a bound). 

The function that takes a Scalar with lower bound I and upper bound u, then 
finds a safe lower bound on I and a safe upper bound on u is implemented as follows: 

enlarge (S 1 u) = S (saf eLowerBound 1) (saf eUpperBound u) 

For convenience we provide a function to convert a number a; to a Scalar with x 
as both lower and upper bound: 

toScalar x = S x x 

7.3. Arithmetic on scalars. We make Scalar an instance of the Mum type class so 
that we can perform arithmetic on scalars (addition (+) , subtraction (-) , negation, 
multiplication (*) and non-negative integer exponentiation (")). 

instance Num Scalar where 

If X € [l,u] then —x G [—u, —I]; negation is exact on safe number so we do not need 
to enlarge the bound: 

negate (S 1 u) = S (-u) (-1) 

If a; € u] then |a;| G [a, b], where a = max{0, 1, —u} and b = — min{0, 1, —u} (it is 
easy to check that this is correct regardless of the signs of I and r). All operations 
involved are exact on safe numbers so we do not need to enlarge the bound: 

abs (S 1 u) = S (maximum xs) (-minimum xs) 
where xs = [0, 1, -u] 

li X g[1,u] and y E [/', u'], then x + y G [I + I' ,u + u']. This operation is not exact 

so we enlarge the bound: 

(S 1 u) + (S 1' u') = enlarge (S (1 + 1') (u + u')) 

If X e [I, u] and y G [/', u'], then x*y € [a,b] where a is the minimum of the numbers 
{1*1' ,l*u' ,u*l' ,u* u'} and b is the maximum of the same numbers. This operation 
is not exact so we enlarge the bound: 

(S 1 u) * (SI' u') = enlcirge (S (minimum xs) (maximim xs)) 
where xs = [1*1 ' , l*u ' , u*l ' , u*u ' ] 

The last two methods are required to complete the implementation of the Num 
instance (f romlnteger provides implicit conversion of integer literals to Scalar 

values): 

fromlnteger = toScalar . fromlnteger 
signum (S 1 u) = error "S.signum: not defined" 
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In order to be able to divide Scalar values using (/) we must also add Scalar 
to the Fractional type class. 

instance Fractional Scalar where 

If X € [l,u] and if /, w have the same sign, then the reciprocal is well defined for x 
and 1/x S I/']- This operation is not exact so we enlarge the bound: 

recip (S 1 u) I l*u > = enlarge (S (1/u) (1/1)) 

I otherwise = error "S. recip: not well-defined" 

The last method is required; it provides implicit conversion of decimal literals to 
Scalar values: 

f romRational = toScalar . f romRational 

7.4. Ordering of scalars. In order to be able to compare Scalar values, e.g. using 
(<), we add Scalar to the Ord type class. If two bounds overlap we declare them 
incomparable and halt the program, otherwise comparison is implemented in the 

obvious way. 

instance Ord Scalar where 

compare (S 1 u) (SI' u') 

I u < 1' = LT 

I 1 > u' = GT 

I 1 == 1' && u == u' = EQ 

I otherwise = error "S. compare: imcompeirable " 

8. Computation with polynomials 

We show how to lift operations on polynomials (of degree d—1) to rectangle sets 
in Mf^ and then how to bound these operations. 

8.1. Representation of polynomials. Polynomials are represented as a list of 
Scalar values (with the first element representing the constant coefficient). Hence 
what we refer to as a 'polynomial' of degree d — 1 is actually a rectangle set in M.'^. 
In this section we lift operations on actual polynomials to such rectangles. We do 
not need to find any bounds on these lifts since this was already done implicitly in 
the previous section. 

8.2. Arithmetic with polynomials. Add polynomials to the Num type class so 
that we can perform arithmetic operations on polynomials. (This implementation 
is a bit more general since it adds [a] to the Num type class for any type a in the 
Mum type class.) 

instance (Num a) => Num [a] where 

Addition: [ci + zqi{z)] + [c2 + zq2{z)] = [ci + C2] + z[qi{z) + q2{z)]. 

(cl:ql) + (c2:q2) = ci + c2 : qi + q2 
[] + p2 = p2 

pi + [] = pi 

Multiplication: [ci + zqi{z)] ■ [c2 + zq2{z)] = [ci ■ C2] + z[ci ■ q2{z) + qi{z) ■ P2iz)], 
where ^2(2) = C2 + zq2{z). 

(cl:ql) * p2@(c2:q2) = cl*c2 : [cl]*q2 + ql*p2 
* _ = [] 
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The remaining methods are straightforward: 

negate p = map negate p 

fromlnteger c = [f romlnteger c] 

abs = error "abs not implemented for polynomials" 

signum = error "signum not defined for polynomials" 

8.3. Polynomial evaluation. Evaluation of the polynomial c + zq{z) at the point 
t is done using the obvious recursion: 

peval (c:q)t=c+t* peval q t 
peval [] _ = 

8.4. Norm of polynomial. We use the £^-norm on polynomials, i.e. ||ao + • • • + 
ttnz'^W = |ao| H h |a„|: 

pnorm = sum . map abs 

8.5. Derivative of polynomial. The derivative of c + zq{z) is implemented using 

the recursion D{c + zq{z)) = q{z) + zDq(z): 

pderiv (c:q) =q+ (0 : pderiv q) 
pderiv [] = [] 

9. Computation with analytic functions 

We show how to lift operations on analytic functions to rectangle subsets in X 
and how to bound these operations. 

9.1. The Function data type. Functions in X are represented as 

f{z)=p{z) + z''h{z), 

where p is a polynomial (not necessarily of degree less than d) and \\h\\ < K, where 
h ^ X. We refer to p as the polynomial part of / and h is called the error of /. 
The value for the degree d is specified in Appendix O 

The Function data type represents an analytic function on the above form (the 
first parameter is the polynomial part, the second parameter is the bound on the 
error): 

data Function = F [[Scalar] ! Scalar deriving (Show,Eq) 

That is. Function represents rectangle subsets of X of the form 

{ao + • • • + a„z" + z'^h{z) \ ak & Ak,k ^ 0, . . . ,n,\\h\\ e /}, 

where {Ak} and / are intervals. Only the upper bound on the error term is needed 
so we do not take care to ensure that the lower bound is correct. Hence, the lower 
bound will be meaningless in general. 

Note that we do allow n > d in the above representation but in general we adjust 
our computations to ensure n < d. We call this operation splitting: if 

f{z) = ao + ••■ + &„ z" + z'^/i(z), 

with n > k > d, then wc can split / at degree k into 

/(z) ao + • • ■ + ak-iz''-^ + z'^iakz''-'^ + ■■■ + anz"'-'^ + h{z)] 

^p'iz) + z''[r{z) + h{z)]. 
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Thus the polynomial part of / after splitting is p' and the error is bounded by 
||r|| + \\h\\ (by the triangle inequality). The implementation of this operation is: 

split k (F p e) = let (p' ,r) = splitAt k p in F p' (e + pnorm r) 

We will now lift operations on analytic functions to the above type of rectangles 
and then find bounds on these operations. 

9.2. Arithmetic with analytic functions. In what follows we let fi{z) = pi{z) + 
z'^hi{z) for i — 1,2, 3, and let /i « /2 = fs where o is the operation under consider- 
ation. 

Make Function an instance of the Num type class so that we can perform arith- 
metic operations on functions (addition (+), subtraction (-), negation, multiplica- 
tion (*) and non-negative integer exponentiation (")). 

instance Num Function where 

Addition of two functions is performed by adding the polynomial part and the error 
separately, — pi+p2 and — h\ + h2, so that ||ft,3|| < ||/ii || + 1|^2|| by the triangle 
inequality: 

(F pi el) + (F p2 e2) = F (pi + p2) (el + e2) 

Multiplication of two analytic functions is given by the equation 

h(z)f2{z) = P1{Z)P2{Z) + Z'' [p^(z)h2(z) + P2{z)hl{z) + z''hl{z)h2{z)\ , 

so that ||/i3|| < |bi||||/i2|| + Ib2||||/ii|| + ||/ii||||/i2||- To ensure that the degree of the 
polynomial part does not increase too much we split it at degree d + Irl 

(F pi el) * (F p2 e2) = split (d+1) (F (pl*p2) e3) 
where e3 = e2*pnorm pi + el*pnorm p2 + el*e2 

The negation of / is —p[z) + z''(— /i(z)) but the error is unchanged since we only 
keep a bound on its norm: 

negate (F p e) = F (negate p) e 

The remaining methods are required to complete the implementation of the Num in- 
stance (f romlnteger provides implicit conversion of integer numerals to Function 
values): 

fromlnteger c = F [f romlnteger c] 

abs = error "abs not implemented for Function" 

signum = error "signum not defined for Function" 

9.3. Norm of analytic functions. The triangle inequality gives ||p(z)+z''/i(z)|| < 
Ibll + \\h\\ (since \z\ < 1): 

norm (F p e) = pnorm p + e 

The norm on the Cartesian product Y — X x X is ||(/,5)11 = ||/|| + H^lj: 
sumnorm (f ,g) = norm f + norm g 

The distance induced by the norm on Y: 
dist (f,g) (f',g') = sumnorm (f-f',g-g') 



We choose to spUt at degree d+1 (instead of the perhaps more natural choice of degree d) 
because the division by z in the definition of the operator T would otherwise cause gs to have 
degree at most d — 2. 
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9.4. Differentiation. The implementation of differentiation of / e X is com- 
plicated by the use of the £^-norm on X, since ||/|| < oo does not imply that 
ll-D/ll < oo. This problem is overcome by only computing the derivative of func- 
tions restricted to a disk of radius strictly smaller than one. That is, we need to 
know a-priori that the function we are differentiating only will be evaluated on 
this smaller disk. Usually we get this information from the fact that we compute 
derivatives like Dfi o /2 and we have bounds on the image of /2 . 

Given / e X we will estimate Df |{|z|<^} where fi < 1. li f{z) = p{z) + z'^h{z), 
then Df{z) = Dp{z) + dh{z)z'^~^ + z'^Dh{z) = pi{z) + z'^hi{z). Here we are faced 
with the problem that we only know the norm of h so all we can say about the 
polynomial part is that pi{z) — Dp{z) 4- sdz'^^^, where s G [— \\h\\]. 

Let h{z) — '^QkZ^, then the error can be crudely approximated as follows: 

fc>i 

Putting all this together we arrive at the following implementation: 

99 deriv mu (F p e) I mu < 1 = F pi el 

I otherwise = error "deriv: mu is not < 1" 
where pi = pderiv p + [S (-s) s] * [0,1] "(d-1) 
el = e / (1 - mu)"2 
s = f romlntegral d * upper e 

Note that fi is passed as a parameter by the caller of this function (it is not a 
constant). As mentioned earlier, usually this function is used to compute expres- 
sions like Dfi o f2 in which case /i will be an upper bound on the radius of a disk 
containing the image of f2- 

9.5. Composition. The implementation of composition of analytic functions /i o 
/2 is split up into two parts. First we consider the special case when /i = pi is a 
polynomial, then we treat the general case. 

Polynomials are defined on all of C so the composition pi o /2 is always defined. 
If pi{z) = c + zq{z) then we may use the recursion suggested by pi o f2{z) ~ 

C + f2{z) -90/2(2:): 

104 compose' (c:q) f2 = (F [c] 0) + f2 * compose' q f2 

The recursion ends when the polynomial is the zero polynomial, in which case 
pi o /2 = 0: 

105 compose ' [] _ = 

In the general case we have to take care to ensure that the image of /2 is contained 
in the domain of /i for the composition to be well-defined. A sufficient condition 
for this to hold is II/2II < 1 since the domain of /i is the unit disk. Under this 
assumption we compute /i o /2 = Pi ° f2 + {f2)'^ • /ii o /2- These two terms are split 
at degree d+lto get pi o f2{z) = pi(z) -I- z'^hi{z) and f2iz)'^ = P2{z) + z''ft,2(2:)0 
Then /i 0/2(2) = p3{z) + z'^hsiz) with ps = pi -l-p2 ■ ^1 0/2 and /13 = ft-i +/12 ■ /ii 0/2. 



See the footnote near the definition of multiphcation of analytic functions for an explanation 
of the choice of degree d+ I. 
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Only the norm oi hi is given so from this we can only draw the conclusion that 
P3{z) ~ pi{z) + s ■ P2iz) for s G [— ||ft.i||, \\hi\\] (s is in fact a function but we may 
think of it as a constant since we are really computing with sets of polynomials). 
The error is approximated using the triangle inequality, Wh^W < \\hi\\ + ||/i2||||/ii||. 

106 compose (F pi el) f2 I norm f 2 < 1 = F (pi ' + [s] *p2 ' ) (el' + el*e2') 

I otherwise = error "compose: |f2| is too large" 
where (F pi ' el') = split (d+1) (compose' pi f2) 
(F p2' e2') = split (d+1) (f2~d) 
s = S (-upper el) (upper el) 

The term s ■ P2{z) can introduce devastating errors into the computation since 
s lies in an interval which has positive upper bound and a negative lower bound 
(if p2 has a coefficient with small error but a large magnitude relative to s, then 
after multiplying with s that coefficient will have an error that is bigger than the 
magnitude of the coefficient). We work around this problem by choosing the de- 
gree d large, since this tends to make the term s smaller. Another way to deal with 
this problem is to include a "general error" term in the representation of analytic 
functions (see [TU]). 

9.6. Derivative of the composition operator. Let S{f,g) — fog, then the 
derivative is given by 

DS(j^g){Sf,Sg)^Dfog.Sg + Sfog. 

Note that when computing Df we must specify as a first parameter the radius of a 
disk strictly contained in the unit disk to which Df is restricted (see Section . 
In the present situation we know that the image of g is contained in a disk with 
radius \\g\\, so Df only needs to be evaluated on the disk of radius ||g||: 

111 dcompose f g df dg = (deriv (norm g) f ~ compose" g) * dg + (df "compose" g) 

9.7. Division by z. If f{z) = aiz + ■ • ■ + a„z" + z'^h{z), then 

f{z)/z = ai + ■ • ■ + a„z"-^ + zJ-^hiO) + z'^h{z), 

where \h{fi)\ < \\h\\ and < \\h\\. Since we do not know the value of h{0) we 
estimate the coefficient it with s G [— 1|^||, ||^||]- We think of this operation as a 
"left shift" , whence the name of this function: 

112 Ishift (F (c:q) e) = F (q + [0,1] '(d-1) * [S (-upper e) (upper e)]) e 
Ishift (F □ e) = F ([0,l]-(d-l) * [S (-upper e) (upper e)]) e 

If the polynomial part of / has a constant coefficient oq ^ then this function will 
not return the correct result, so we take care to only use it when we know that 
ao = 0. 

9.8. Point evaluation. If f{z) = p{z) + z'^h{z), then f{t) = p(t) +t'^-s for some 
s G [— \\h\\]. We also check that t is in the unit disk otherwise the program is 
terminated with an error: 

114 eval (F p e) t I abs t < 1 = peval p t + t~d * (S (-upper e) (upper e)) 
I otherwise = error ("eval: not in domain t=" ++ show t) 

Note that the further away t is from 0, the more error is introduced in the evaluation. 
For t = the error term has no influence on the evaluation. 
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9.9. Scaling. As a convenience we define operators to scale an analytic function 
by a scalar on the on the right. The precedence for these operators are the same 
as for their 'normal' counterparts. 

Multiplication satisfies [p{z) + z'^h{z)] ■ x — x ■ p{z) + z'^[x ■ h{z)] and division is 
handled similarly. Note that the error term is affected: 

116 infixl 7 .*, ./ 

(F p e) .* X = F (p * [x]) (e * abs x) 
(F p e) ./ X = F (p * [1/x]) (e / abs x) 

9.10. Polynomial approximation. Let f{z) — p{z) + z'^h{z). To approximate 
/ by a polynomial we first discard the error term z'^h{z), then we disregard the 
errors in the coefficients of p. That is, for p{z) = aq + • ■ • + o„2:" with Uk € [a^, a^] 
we replace with the mean hk = (a^T + )/2 (we 'collapse' the bounds on a^). 
Finally, we lift this operation to pairs of functions: 

119 approx (f,g) = (approx' f, approx ' g) 

where approx' (F p _) = F (map (toScalar . collapse) p) 
collapse (S 1 u) = (l+u)/2 

9.11. Operator norm. Let £,k{z) — z^ so that {Cfe}fc>o is a basis for X. A basis 
for Y is {?/fc}/c>o, where ■q2k = (Cfe,0) and ?72fe+i = (0,^fe). This set is implemented 
as follows: 

122 basis = interleave (zip basis' (repeat 0)) (zip (repeat 0) basis') 
where basis' = map xi [0..] 

xi k = F (replicate k ++ [1]) 

Proposition 9.1. If L :Y is a linear and bounded operator, then 
\\L\\ = max{||L7yo||, . . • , ||Lr,2d-i||, sup \\L{h, 0)||, sup ||L(0, 

heBa h£Ba 

where Bd {z'^h{z) \ \\h\\ < 1}. 

This is a consequence of using the .^^-norm on X. 

Given a linear operator op acting on a list of tangent vector^, we estimate the 
operator norm by applying it to the first 2d basis vectors and to the sets Bd x 
and X Bd- Then we compute the upper bound of the norm of the results and take 
the maximum: 

125 opnorm op = maximum $ map (upper . sumnorm) $ op tangents 

where tangents = (F [] 1,0) : (0,F [] 1) : take (2*d) basis 

Note that Bd is represented by the set of functions with no polynomial part and an 
error bounded by 1, which is the same as F [] 1. 



The linear operator acts on a sequence of tangent vectors since this is how we have imple- 
mented the derivative of the main operator. 
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9.12. Construction of balls. We cannot exactly represent arbitrary balls in X 
with the Function type. Instead we construct a rectangle set which is guaranteed 
to contain the ball. 

Thus, a bound on a ball of radius r centered on an analytic function (in our case 
it is always a polynomial, i.e. e=0) can be implemented as follows: 

127 ball r (f ,g) = (ball' r f, ball' r g) 

where ball' r (F p e) = F (map (+ S (-r) r) p) (e + toScalar r) 

9.13. Newton's method. This is our variant of Newton's method on Y: 

{f,g)^iM-I)-\M-T){f,g), 

where M is a, 2d x 2d matrix passed as the first parameter. The second parameter 
is T{f,g) and the third parameter is (/,(?). When lifting M into Y we project the 

error term to zero by letting s = and when lifting (M — we preserve the 
error term by letting s = 1 (see below for how this lifting is done): 

129 newton m (tf ,tg) fg = fg' 

where (mf ,mg) = liftPolyOp (apply m) fg 

fg' = liftPolyOp 1 (solve $ subtractDiag m 1) (mf -tf ,mg-tg) 

Let / = p{z) + z'^h{z) where degp < d and let A be the linear operator repre- 
sented (in the basis {z'^}) by the infinite matrix 

fM 
\0 si 

where M is -a d x d matrix and I is the infinite identity matrix. We lift the linear 
operator A into X by 

Af{z) = Mp{z) + z'^{s-h{z)). 

The following function implements this lifting into Y. We split / and g to ensure 
that their degrees are at most d — 1, and since our linear algebra routines require 
its input in one vector we interleave the polynomial parts. Also, instead of passing 
M we pass a linear operator op which allows us to use one function to lift matrix 
multiplication (apply) and solution of linear equations (solve): 

132 liftPolyOp s op (f,g) = (F pf (s*ef), F pg' (s*eg)) 

where fg@(F pf ef, F pg eg) = (split d f , split d g) 

(pf',pg') = uninterleave $ op (interleavePoly fg) 

When interleaving the polynomial parts of two functions we first pad the polyno- 
mials with zeros to ensure their lengths are exactly d (e.g. oq + aiz is padded to 
flo + fliz + Oz^ + • • ■ + Oz'^^^). Hence the resulting vector always has length 2d: 

135 interleavePoly (F p _, F q _) = interleave (pad p) (pad q) 
where pad x = tsike d $ x ++ (repeat 0) 

Appendix A. Linear algebra routines 

In this section we implement a simple linear algebra library to compute matrix- 
vector products and to solve linear equations. 

A matrix is represented as a list of its rows and a row is a list of its elements. 
A vector is just a list of elements (we think of them as column vectors). This is 
a very simplistic library so no checking is done to ensure that matrices have the 
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correct dimensions (e.g. it is quite possible to create a 'matrix' with rows of differing 
lengths). 

A.l. Matrix-vector product. Computing the matrix-vector product Mx is fairly 
straightforward: 

137 apply m X = map (dotProduct x) m 

The dot product of vectors a and b: 

138 dotProduct a b = sum $ zipWith (*) a b 

Note that if a and b have different lengths, then the above function will treat the 
longer vector as if it had the same length as the shorter. 

A. 2. Linear equation solver. The following function solves the linear system of 
equations Mx = b. It is a simple wrapper around a function which solves a linear 
system given an augmented matrix. 

139 solve m b = solveAugmented $ augmentedMatrix m b 

The augmented matrix for M and b is [M 6] , i.e. the matrix with b appended as 

the last column of M: 

140 augmentedMatrix = zipWith (\x y -> x ++ [y] ) 

We now implement the linear equation solver which takes an augmented matrix 
as its only parameter. It is implemented using Gaussian elimination with partial 
pivoting. The only novelty compared with a traditional imperative implementation 
is that we solve the equations recursively. 

Given a n x (n + 1) augmented matrix M' first perform partial pivoting, i.e. the 
row whose first element has the largest magnitude is moved to the top to form 
the matrix M. Assuming that we already have the solution for 0:2, ■ . • , a;„ we can 
compute Xi = (mi_„+i — X]j=2 ^^'^ done. The solution for 

X2, - ■ ■ ,Xn is found recursively as follows: perform a Gaussian elimination on M to 
ensure that all rows except the first start with a zero to get a matrix M. Throw 
away the first row and column of M to get a (n — 1) x n matrix N' and solve the 
linear system with augmented matrix A'"'. The solution to this system is X2^ * • • ? ^ri' 

141 solveAugmented [] = [] 

solveAugmented m' = (last mlt - dotProduct mlt x) / mil : x 
where mS((mll:mlt) :_) = partialPivot m' 
X = solveAugmented $ eliminate m 

Partial pivoting is done by first finding a list of all possible ways to split the ma- 
trix M into a top and a bottom half. This list is searched for the split which has a 

maximal first clement in the bottom half. The maximal split is then reassembled 
into one matrix by moving the top row of the bottom half to the top of the matrix. 

145 partialPivot m = piv:mtop ++ mbot 

where (mtop,piv:mbot) = meiximiimBy comparePivotElt (splits m) 

comparePivotElt (_,(a:_):_) (_,(b:_):_) = compaire (abs a) (abs b) 

The following routine uses Gaussian elimination to ensure that the first element of 
all rows except the first starts with a zero. That is, we add a suitable multiple of 
the first row to the other rows one at a time: 
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148 eliminate ( (mil :mlt) :mbot) = foldl appendScaledRow [] mbot 

where appendScaledRow a (r:rs) = a ++ [scaleAndAdd (-r/mll) mlt rs] 
scaleAndAdd s a b = zipWith (+) (map (*s) a) b 

Remark A.l. When using the above hnear equation solver with matrices and vectors 
over intervals (of type Scalar) there is a question of what the 'solution' represents. 
As always, we are computing bounds on solutions: if M is in some rectangle set 
[M] of matrices and b is in some rectangle set [b] of vectors, then the above routine 
will compute a rectangle set [x] such that if x is a solution to Mx = b then x £ [x]. 

Note that our solver will compute rather loose bounds on the solution set, see 
e.g. [S] for ways of finding sharper bounds. 

Appendix B. Supporting functions 

Given a square matrix M and a number x compute M — xl, i.e. subtract x from 
every diagonal element of M: 

151 subtractDiag m x = foldl f [] (zip m [0..]) 

where f m' (r,k) = let (h,t:ts) = splitAt k r 
in m' ++ [h ++ [t-x] ++ ts] 

Given a list, return all possible ways to split the list in two: 

154 splits X = splits' [] X 

where splits' _ [] = [] 

splits' X y@(yh:yt) = (x,y) : splits' (x ++ [yh] ) yt 

Interleave two lists a and b, i.e. construct a new list by taking the first element 
from a, then the first element from b and repeating for the remaining elements. 

157 interleave a b = concat $ zipWith (\x y -> [x,y]) a b 

Perform the 'inverse' of the above function, i.e. take a list c and construct a pair 
of lists (a,b) such that interleave a b = c: 

158 uninterleave = unzip . pairs 

Given a list, partition it into pairs of adjacent elements: 

159 pairs [] = [] 

pairs (x:y:rest) = (x,y) : pairs rest 

pairs _ = error "list must have even length" 

Appendix C. Input to the main program 

The degree of the error term in our representation of analytic functions: 

162 d = 13 : : Int 

The radius of the ball on which $ is a contraction: 

163 beta = 1 . Oe-7 : : Double 

The radii for the domains of and ip: 

164 sf = 2 . 2 : : Scalar 
sg = . 5 : : Scalar 

The initial guess for the fixed point: 

166 guess = (F [-0.75, -2.5] 0, F [6.2,-2.1] 0) 
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Appendix D. Running the main program 

This document contains all the Haskell source code needed to compile the pro- 
gram into an executable. Given a copy of the I^Tg^X source of this document 
(assuming the file is named Imca. Ihs), use the following command to compile itF^ 

ghc — make -02 Imca. Ihs 

This produces an executable called Imca (or Imca.exe if you are using Windows) 
which when called will execute the main function. 
Here is the output of running the main program: 

radius = l.Oe-7, 

|Phi(f)-f I < 4.830032057738462e-9, 

iDPhil < 0.1584543808202988 

This output was taken from a sample run using GHC 6.12.1 on Mac OS X 10.6.2. 
The nmning time on an 1.8 GHz Intel Core 2 Duo was less than 10 seconds. 

Appendix E. Haskell mini-reference 

This section introduces some of the features and syntax of Haskell to help any- 
body unfamiliar with the language read the source code. It is assumed that the 
reader has some prior experience with an imperative language (Java, C, etc.) but 
is new to functional programming languages. Table [1] below collects examples of 
Haskell syntax used in the source code and can be used to look up unfamiliar expres- 
sions. For more information on the Haskell language go to http : / /haskell . org. 

Haskell is a, functional language. Such languages differ from imperative languages 
in several significant ways, e.g.: there are no control structures such as for loops 
and data is immutable so there is no concept of variables (memory locations) that 
can be written to. 

Basic types include: booleans (True, False), numbers (e.g. -1, 2.3e3, integers 
of any magnitude are supported), tuples (e.g. (1, 'a' ,0.3), elements can have 
different types), and lists (e.g. [1,2,3], all elements must have the same type). 
Functions are on the same level as basic types so they can e.g. be passed as param- 
eters to other functions. 

Functions are defined like f parameters = expression where f is the function name 
and there can be zero or more parameters. Note that there are no parentheses 
around parameters and that parameters are separated by spaces. Function calls 
have very high precedence, so f x~2 is the same as (f x)~2, not f (x~2). The 
keywords let . . in and where can be used to bind expressions to function-local 
definitions (i.e. local functions or variables). 

New data types can be defined using the data construct. For example, data 
Interval = I Double Double defines a type called Interval which consists of two 
double precision fioating point numbers (i.e. the endpoints of the interval). New 
values of this type are constructed using the value constructor which we called 1, 
e.g. 10 1 defines the unit interval. 

Functions can be defined with pattern matching on built-in and custom data 
types. For example, len (1 a b) = abs (b-a) defines a function len which re- 
turns the length of an interval (for the custom data type Interval). 



The GHC compiler can be downloaded for free from http://haskell.org. 
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We often use pattern matching on lists, where [] matches the empty hst and 
(x:xs) matches a hst with a least one element and binds the first element to x and 
the rest to xs (read as plural of x). The notation v@(x:xs) can be used to bind 
the entire list to v on a match. 

The notation _ may be used to match anything without binding the match to 
a variable, e.g. firstZero (x:_) = x == defines a function which returns True 
if the first element of a non-empty list is equal to zero (and throws an exception if 
called on the empty list [] ). 

Type classes are a way of declaring that a custom data type supports a certain 
predefined collection of functions and also allow for 'overloading' of functions (and 
operators, which can be turned into functions as noted in the example for (+) 
in Table [T]). We only mention type classes because we come across them when 
implementing Scalar and Function. The pre-defined type classes we use are Num 
(for (+), (-), (*), (~), abs). Fractional (for (/), (~~)), Eq (test for equality), 
and Show (for conversion to strings). 



Table 1. Examples of Haskell syntax used in the source code. 



Expression 


Description 


fl X = 2*x 


define a function f 1 which doubles its argument 


f2 X y = x+y 


define a function f 2 which adds its two arguments 


fl 3 


apply f 1 to 3 (=6) 


f2 3 4 


apply f 2 to 3 and 4 (= 7) 


f2 2 (fl 3) 


apply f 2 to 2 and 6 (the result of f 1 3) (= 8) 


f2 2 $ fl 3 


same as above (the operator $ is often used in this way to 




avoid overuse of parentheses) 


f2 2 fl 3 


error (this means: compute f2 2 f 1 and apply the result 




to 3, but 2+f 1 does not make sense) 


\x -> 2*x 


define the anonymous function x ^ 2x 


f2 3 


apply f 2 to 3 (= the function \x -> 3+x) 


fl . f2 3 


function composition (= the function \x -> 2* (3+x)) 


3 ~f2~ 4 


turn function (in backticks) into an operator (—7) 


(+) 3 4 


turn operator (in parentheses) into a function (=7) 


(3*) 


fix first parameter to 3 (= the function \x -> 3*x) 


g X 1 x<0 = -1 


define the sign function g using guards (the I symbols) 


1 x>0 = 1 




1 x==0 = 




num [] = 


define a function which counts the number of elements in 


num (_:xs) 


a list using pattern matching 


= 1+num xs 




[1,2] 


a list (all elements must have the same type) 


[] 


the empty list 


[1. .] 


the list of all positive integers 


[2. .5] 


list enumeration with bounds (=[2,3,4,5]) 


1 : [2,3] 


append element to beginning of list (=[1,2,3]) 



Continued on the next page. . . 
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Expression 



Description 



[1,2,3] ! ! 
[1,2] ++ [3] 
'a' 

"abc" 
2-3 

2"(-l) 
Ca' ,2) 
fst Ca' ,2) 
snd Ca' ,2) 
map f 1 [1. .] 
[f2 a b I 

a <- [1,2] , 
b <- [3. .5]] 
foldl ±2 1 [3,5] 
iterate fl 1 
maximum [1,4,2] 
maximumBy f x 
minimum [1,4,2] 
splitAt 2 [1,4,2] 
take 3 [7. .] 
zip [1..] [3,4] 
unzip [(1,3), 
(2,4)] 
zipWith f2 [1. .] 
[3,4] 

repeat 
replicate 3 
sum [3,-1,4] 
transpose m 
putStrLn "hi" 
show 1.2 
error "ohno" 



access list elements by zero-based index (=1) 
concatenate two lists (=[1,2,3]) 
a character 

a string, i.e. a list of characters (=['a' , 'b' , 'c']) 

non-negative integer exponentiation (=8) 

integer exponentiation (=0 . 5) 

a pair (the elements need not have the same type) 

access first element in a pair (= ' a ' ) 

access second element in a pair (=2) 

apply f 1 to all elements in the list (=[2,4,8, . .]) 

list comprehension, i.e. {/2(a, 6) | a G {1,2}, 6 G {3,4,5}} 

(= [4,5,6,5,6,7]) 

fold left over list (compute f 2 1 3 = 4, then f 2 4 5) (=9) 

compute orbit of 1 under fl (=[1,2,4,8, . .]) 

return maximum clement in a list (=4) 

as above, but using f to compare elements of the list x 

return minimum element in a list (=1) 

split list in two at given index (=( [1,4] , [2] )) 

take the first 3 elements from the list (=[7,8,9]) 

join two lists into a list of pairs (=[(1,3), (2, 4)]) 

'inverse' of zip (= ( [1 , 2] , [3 , 4] ) ) 

like zip, but use f 2 to join elements (= [4,6] ) 

infinite list with one element repeated (=[0,0,0, . . ] ) 
finite list with one element repeated (=[0,0,0]) 

sum elements in list (=6) 

the transpose of the matrix m (m is a list of lists) 
print hi to standard out and and append a new line 
turn the number 1 . 2 into the string "1.2" 
abort program with error message ohno 
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